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The high-resolution wave-propagation method for computing the nonhydrostatic 
atmospheric flows on meso- and micro-scales is described. The design and implementation 
of the Riemann solver used for computing the Godunov fluxes is discussed in detail. The 
method uses a flux-based wave decomposition in which the flux differences are written 
directly as the linear combination of the right eigenvectors of the hyperbolic system. The 
two advantages of the technique are: 1) the need for an explicit definition of the Roe matrix 
is eliminated and, 2) the inclusion of source term due to gravity does not result in 
discretization errors. The resulting flow solver is conservative and able to resolve regions of 
large gradients without introducing dispersion errors. The methodology is validated against 
exact analytical solutions and benchmark cases for non-hydrostatic atmospheric flows. 


I. Introduction 

T he high-resolution wave-propagation method was introduced by LeVeque as a general algorithm for solving 
hyperbolic conservation laws (LeVeque 1996; LeVeque 1997; LeVeque 2002). It has been successfully applied 
to a wide range of physical problems (LeVeque 2002). Traditionally, atmospheric modelers have favored the 
centered finite difference schemes. These types of schemes exhibit non-physical spurious oscillations, which can 
contaminate the numerical results and introduce false negatives in microphysical scalars. At smaller spatial and 
temporal scales, large gradients of velocities and other physical quantities can develop, which create stability 
problems for such schemes. The high-resolution wave-propagation methods are conservative finite volume 
discretizations and have the ability to resolve regions of steep gradients accurately, thus avoiding dispersion errors in 
the solution. Positivity of scalars (an important factor when considering the transport of microphysical quantities) is 
also guaranteed by applying the total variation diminishing (TVD) condition appropriately. 

The conservation laws used in atmospheric models are simplified by assuming that the atmospheric processes are 
adiabatic - the conservation of total energy is replaced by the conservation of potential temperature (Das 1979). 
Ooyama (1990) stated that ’’The first law of thermodynamics is a statement of conservation in the processes that 
involve more forms of energy than the purely mechanical. Since atmospheric processes are predominantly 
adiabatic, it is appropriate to express the first law in terms of entropy conservation.” This has often been the 
argument for using the potential temperature conservation equation rather than the total energy. Some examples of 
the flow model with the potential temperature conservation can be found in Proctor (1987) and Bacon et al. (2000). 
The operational mesoscale models such as the Regional Atmospheric Modeling System (RAMS) developed by 
Pielke et al. (1993) and the Coupled Ocean/ Atmosphere Mesoscale Prediction System (COAMPS) developed by 
Hodur (1997) are all based on equations which express the conservation of total energy in terms of functions for 
potential temperature. The conservation laws with total energy equation are rarely used for atmospheric flow 
computations. Some exceptions include, Botta et al. (2004), Ahmad et al. (2007) and Giraldo and Restelli (2008). 

In this study, the high-resolution wave-propagation method is applied to meso- and micro-scale atmospheric 
flows. The equation set comprises of conservation of total energy rather than potential temperature. There are two 
major advantages (Giraldo and Restelli 2008) in using the conservation of total energy equation: 1) energy is 

conserved to machine precision and, 2) true viscous stresses are accounted for if the diffusion term is included. In 
the following sections the governing equations and the numerical technique are described. The results of benchmark 
simulations performed to evaluate the accuracy of the numerical scheme and its applicability to nonhydrostatic 
atmospheric flows are discussed in detail. 
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II. Governing Equations 

The basic equations governing fluid flow comprise of a set of partial differential equations for the conservation 
of mass, momentum and energy. In two dimensions they can be written as: 
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p is the density of fluid, u is the velocity component in the x-direction, w is the velocity component in the z-direction 
and p is the pressure. 

E is the total energy per unit volume: 


E = pe + /? — (w 2 + w 2 ) 


(3) 


where, — {u 2 + w 2 ) is the specific kinetic energy. The specific internal energy, e is given by an equation of state: 


(/-!)/> 


(4) 


In the above relations, y is the ratio of specific heats and g is the acceleration due to gravity. In this study the 
atmosphere is assumed to be dry and the source term S is set to zero (S would represent heat sinks and sources due to 
microphysical processes.). The viscous fluxes, F y and H y are given by: 
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In Eq. (5) p is the dynamic viscosity, T is temperature, c p is the specific heat at constant pressure and Pr is the 
Prandtl number. 

The viscous shear stress tensor, Zjj is given by: 
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III. Numerical Scheme 

The conservation laws in Eq. (l)-(2) without the source terms can be written in the discrete form in ID as: 
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where U is the vector of conserved quantities, and F is the vector of inter-cell fluxes calculated at the control surface 
of each control volume using either an exact or approximate Riemann solver. At and Ax are the time step and mesh 
resolution in x-direction, respectively. LeVeque (2002) and Bale et al. (2002) suggest using a flux-based wave 
decomposition, in which the flux differences F. (U i ) - F._, (f/ M ) are written directly as a linear combination of the 

right eigenvectors, r ?_ x/2 : 
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A-i/2=^r-i/2(^(V)-^-i(V-i)-Ax^_ 1/2 ). 
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The vectors = J3 p r p are called /-waves and contain flux increments rather than increments in U. R i _ x/2 is the 
matrix of right eigenvectors. \j/ - pg , is the source term due to gravity. Eq. (7) can now be re-written as: 
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The fluctuations H A U i _ l/2 and H + Af/ M/2 contribute to the cell-averaged quantity f// due to the wave 

propagation across the cell interfaces. In the above relations, /^ 1/2 are the wave speeds given by the eigenvalues 

of the hyperbolic equation set. Higher-order accuracy in space can be achieved by adding a correction term 
(LeVeque 1996: LeVeque 1997; LeVeque 2002): 
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and, Z p is the limited value of Z p . Given the /-waves and the wave speeds, the flux differences can be computed 
by summing up the left and right going waves across a cell interface. In the above relations, the sweep in x-direction 
is implied. Similar methodology can be used for computations in the z-direction with the addition of a source term 
in the z-direction. The quantities on cell faces are calculated by taking the average of cell-centered quantities on the 
either side of the face. 


IV. Results 

In this section the numerical scheme is evaluated against several benchmark cases - (1) isentropic vortex 
evolution (Yee at al. 1999); (2) transport of density distribution (Deiterding 2000); (3) Rayleigh-Taylor instability 
(Almgren 2010); (4) inertia-gravity waves on non-hydrostatic scale (Skamarock and Klemp 1994); and (5) density 
current (Straka et al. 1993). In the first four cases the Euler equations are computed and in the last test case the 
Navier-Stokes equations are computed. The first two test cases have exact analytical solutions. 

A. Isentropic Vortex 

The evolution of a stationary and a moving inviscid isentropic vortex (Yee et al. 1999) is described in this 
section. In this case, there is no source term due to gravity. The computational domain was defined as 
(x,z) e [0,10] x [—5,5] m with t e [0,100] s. The mesh size was set to 100 x 100 cells. Periodic boundary conditions 

were used in all directions for the moving vortex and farfield conditions were used for the stationary vortex. For the 
stationary vortex, the u-ve locity and w-velocity were set to zero and in the case of the moving vortex the i/-ve locity 
was set to 0.2 m/s. The initial density and pressure were set to Ikg/m 3 and 1 Pa respectively. Perturbations in the 
flow field (Su,Sw) and density ( Sp ) were added such that SS = 0 ? where, S is the entropy (Yee at al. 1999): 

(Su,Sw) = — e (1 ~ r )/2 (-z,x) (15) 


The initial density field was defined by: 
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In the above equations (5 is the vortex strength which was set to 5 . r 2 = x 2 + z 2 and (x, z) = (x - x 0 , z - z 0 ) . 
The initial vortex location is given by (x 0 , z 0 ) = (5,0) . The pressure was initialized as follows: 

P = P r (17) 

The root mean square errors in density (£p) for the stationary vortex case were calculated for different mesh 
resolutions and are listed in Table 1. The computed density fields and the comparison with exact solutions along the 
mesh centerline are shown in Figure 1 (stationary vortex) and Figure 2 (moving vortex), respectively. The 
computed solutions are in good agreement with the exact solution and exhibit little numerical diffusion. The 
solution on a higher-resolution mesh (200 x 200 cells) and a comparison with the exact solution along the mesh 
centerline is shown in Figure 3. 


Table 1: Stationary Isentropic Vortex 


# of cells 

Av(m) 

E P 

50x50 

0.2 

9.41E-03 

100 X 100 

0.1 

1.34E-03 

200 x 200 

0.05 

1.82E-04 
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X (m) x M 


Figure 1. Stationary Vortex. Computed density (k g/m 3 ) at time =100s (left) and the comparison of the 
computed solution with the analytical solution at z = 0 and time = 100s (right). Mesh size = 100 x 100 cells. 
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Figure 2. Moving Vortex. Computed density (kg/m 3 ) at time =100s (left) and the comparison of the computed 
solution with the analytical solution at z = 0 and time = 100s (right). Mesh size = 100 x 100 cells. 
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Figure 3. Moving Vortex. Computed density (kg/m 3 ) at time =100s (left) and the comparison of the computed 
solution with the analytical solution at z = 0 and time = 100s (right). Mesh size = 200 x 200 cells. 


5 

American Institute of Aeronautics and Astronautics 





B. Transport of Smooth Density Distribution 

Deiterding (2000) suggested a simple test case with an analytical solution for evaluating numerical schemes for 
the Euler equations with a gravitational source term. The Euler equations can be reduced to a single transport 
equation for density: 


u(x,z,t) = u 0 
w(x,z,t) = w 0 -gt 
p(x,z,t) = p 0 


The following relation describes the advection of an initial arbitrary density profile (Deiterding 2000): 


p(x,z,t) = p 0 


1 2 

x-u 0 t , z -w 0 t + — gt 


(18) 
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The computational domain in this case was defined as (x,z) e [0,2] x [0,2] m with te [0,0.5] s. Periodic 
boundary conditions were used in all directions. The transport of density distribution is given by: 


x c (t ) = x 0 +u Q t 

1 2 

z c (t) = z 0 +W 0 t--gt 


The initial density distribution is defined as follows: 
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where, 


r(x, z, t) = ^j(x-x c (t)f + (z-z c (t)f (22) 

In the above equations p\ = Ikg/m 3 , pi - 0.05 kg/m 2 ’, x 0 = 0.75m, z 0 = 0.75 m, u 0 = 1 m/s, w 0 - 1.25 m/s,po = IPa, g 
= hn/s 2 , and y= 1.39. 

The simulation was run for different mesh resolutions and the root mean square errors in the density (E p ) and the 
velocity (E Y ) fields were calculated for each simulation (Table 2). 

The density distribution at time = 0.5s is shown in Figure 4. The computed density profile along the mesh 
diagonal is also compared with the exact analytical solution in Figure 4. The computed results are in excellent 
agreement with the exact analytical solution. 


Table 2: Transport of Smooth Density Distribution 


# of cells 

Av(m) 

E P 

E y 

40x40 

0.05 

0.536237E-03 

0.289277E-15 

80x80 

0.025 

0.210343E-03 

0.361827E-15 

160 x 160 

0.0125 

0.940737E-04 

0.362340E-15 

320 x 320 

0.00625 

0.4490 19E-04 

0.776312E-15 
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Figure 4. Transport of Density Distribution. The computed density distribution ( kg/m 3 ) at time = 0.5s is 
shown in the left panel and the comparison of the computed solution with the analytical solution along the 
mesh diagonal at time = 0.5s is shown in the right panel. 




C. Rayleigh-Taylor Instability 

Simulation of the Rayleigh-Taylor instability has been described in several publications (Almgren et al. 2010; 
Liska and Wendroff 2003; Li et al. 1996) and is often used as a validation test case. The computational domain in 
this case was defined by (x,z) e [0,0.5] x [0,1] m with t e [0,2.5] s. The mesh size was 256 x 512 cells. Periodic 
boundary conditions were used in the lateral and the top and bottom boundaries were set to solid walls. The initial 
w-velocity and the w- velocity were set to zero. The density was set to pi = \kg/m 3 in the lower half of the domain 
and P 2 = 2 kg/m* in the upper half of the domain. Pressure was initialized using the hydrostatic equation: 


p(z) = 


| Po ~ P\S Z 

1^0 - P\g L z I2-Pig(z-L z l 2) 


z<LJ 2 
z>LJ2 


(23) 


where p 0 is the reference pressure and was set to 5 Pa. The acceleration due to gravity, g was set to 1 mis 2 . L z = Im 
is the height of the computational domain. A single-mode perturbation was introduced in density at the interface of 
heavier and lighter fluids: 


P(x,z) = p l 


_l_ Pi P\ 
2 


1 + tanh 


0.005 J 


where, 


¥(x) = L + 0,01 cos (4 *) + cos(4 n (L x -x)) 
2 2 


(24) 


(25) 


In Eq. (25), Z x = 0.5m is the width of the computational domain. The tanh function in Eq. (24) is used to 
smooth the density profile at the interface (Almgren et al. 2010). The initial conditions and the evolution of the 
instability are shown in Figure 5. The scheme was able to capture the onset of the instability accurately and the 
solution is in good agreement with previously reported results in the literature (e.g., Almgren et al. 2010). 
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Figure 5. Rayleigh-Taylor Instability. The initial density field (top left panel); density field at Is (top right 
panel); at 2s (bottom left panel); and at 2.5 s (bottom right panel). Mesh size = 256 x 512 cells. 
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D. Non-hydrostatic Inertia-Gravity Waves 

The simulation of inertia-gravity waves (Skamarock and Klemp 1994) on the non-hydrostatic scale is described 
in this section. The computational domain was defined by (x, z) e [0,300] x [0,10] km with t e [0,3000] s. The mesh 

had a resolution of 1 km in the x and 50 m in the z-direction (300 x 200 cells). Periodic boundary conditions were 
used in the lateral and the top and bottom boundaries were set to solid walls. The domain was initialized by a 
constant Brunt-Vaisala frequency, N = 10'V 1 . The waves were excited by an initial potential temperature 
perturbation given by: 


9(x,z, 0) = A9 0 


sin [nz! H) 

1 + (x -x c ) 2 la 2 


(26) 


The amplitude of the initial potential temperature perturbation, A 6^ was set to 10 ~ 2 K. The height H of the 

domain was 10 km, the perturbation half width was a = 5 km. The perturbation was initialized at x c = LJ 3, where, L x 
is the width of the domain (300 km). A uniform ^-velocity of 20 m/s was imposed in the domain and the w-velocity 
was set to zero. The initial potential temperature perturbation field is shown in Figure 6 and Figure 7 shows the 
computed potential temperature perturbation at time = 3000s. Overall, the solution agrees well with other results 
reported in literature (Ahmad and Lindeman 2007; Giraldo and Restelli 2008). Compared to Giraldo and Restelli 
(2008) who used a higher-order Discontinuous Galerkin Method, the results are slightly diffusive (the wave in the 
middle of the domain has a lower peak). The final minimum and maximum potential temperature perturbation in 

Giraldo and Restelli (2008) are O' e [— 1.51xl0 -3 , 2.78 xlO -3 ] K compared to O' e [-1.41x1 0 -3 , 2.83 xlO -3 ] K m the 
current study. 
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Figure 6. Inertia Gravity Waves. The initial potential perturbation field ( K) is shown in the figure. The field 
minimum = 0 K and the maximum = 0.01IT. Mesh size = 300 x 200 cells. 



Figure 7. Inertia Gravity Waves. The potential temperature perturbation field (K) at time = 3000s. The field 
minimum = -0.001417T and the maximum = 0.002837T. Mesh size = 300 x 200 cells. 
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E. Density Current 

The simulation of density current (Straka et al. 1993; Ahmad and Lindeman 2007; Giraldo and Restelli 2008) is 
described in this section. The computational domain was defined by (x, z) e [0,25] x [0,6.4] km with t e [0,900] s. 

The mesh had a uniform resolution of 50m in both x and z-directions (500 x 128 cells). Farfield/outflow boundary 
conditions were used for the right lateral boundary and all other boundaries were set to solid walls. The domain was 
initialized for a neutral atmosphere by setting the potential temperature at 300 K. The initial i/-velocity and the w- 
velocity were set to zero. A vertical temperature profile was defined for the entire domain using the following 
relation: 


T = 300- 


zg 


(27) 


A cold bubble was initialized by adding a perturbation in the temperature field using the following relation: 


AT = < 


0.0 

-15.0 


cos(/rZ)+1.0 


if L >1.0, 
if L <1.0 


(28) 


where, 



where, (x c , z c ) = (0,3 ) km and (x r ,z r ) = (4,2) km. The initial potential temperature perturbation field is shown in 

Figure 8. A constant viscosity (// = 75 m 2 /s) was used in the calculation of momentum and energy viscous fluxes. 
The Prandtl number was set to 0.76. 

This test case can be considered as an idealized downburst gust front. As the cold air descends due to negative 
buoyancy, strong downdrafts develop at the center of the cold bubble. When the cold air reaches the ground, it is 
rolled up and a front is formed. As this front propagates, shear is generated at the top boundary of the front. The 
benchmark solution (Straka et al. 1993) consists of three rotors, which develop at the top boundary of the front due 
to Kelvin-Helmholtz type instability. The potential temperature perturbation field at time = 900s is given in Figure 
9, which shows the three characteristic rotors of the Straka benchmark. The front is located at x = 15326m. The 
location of the front in the Straka reference solution is at x = 15537m. The location of front in Giraldo and Restelli 
(2008) is at x = 14789m. The reference solution in Straka et al. (1993) and Giraldo and Restelli (2008) are on a 25m 
resolution mesh whereas a mesh resolution of 50 m was used in the current study. The present results also compare 
well with other reported solutions of the benchmark (e.g., Janjic et al. 2001; Xue et al. 2000). 



x(m) 


Figure 8. Density Current. The initial potential perturbation field ( K) is shown in the figure. The field 
minimum = -16.6 K and the maximum = OK. Mesh size = 500 x 128 cells. 
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x(m) 

Figure 9. Density Current. The potential temperature perturbation field (A) at time = 900s. The contour 
values shown in the figure are from -9 K to 0 K. The Mesh size = 500 x 128 cells. 

V. Summary 

LeVeque’s high-resolution wave-propagation method was applied to non-hydrostatic atmospheric flows on 
meso- and micro-scales. The equation for total energy conservation was used instead of potential temperature 
conservation in the simulations. The results show that the numerical method is fully conservative and is able to 
discretize the source due to gravity without introducing errors in the solution. Comparisons with exact analytical 
solutions show that the flow solver has relatively small diffusion and dispersion errors. In the case of isentropic 
vortex and density transport case, the computed solutions are in excellent agreement with the exact analytical 
solutions. The simulations of the Rayleigh-Taylor instability, nonhydrostatic inertia gravity waves, and the density 
current also compare well with previously reported results in the literature. Since the eigenstructure of the Euler 
equations is well known, the scheme can be extended to three dimensions for simulating realistic flows. It is 
especially well-suited for simulating atmospheric flows in which sharp gradients of temperature, potential 
temperature and winds can develop (e.g., fronts, supercell thunderstorms, tornadoes). 
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